1. LECTURA DE DADES

En aquest document farem el mateix estudi que amb les dades de Mallorca però utilitzant les d’Eivissa i Formentera. Els càlculs i procediments són anàlegs al de l’illa de Mallorca, per això obviarem alguns detalls.

En primer lloc, grafiquem en forma de sèrie temporal aquestes dades:

turisme <- read_excel("IBESTAT.xls")
turisme$Data <- gsub("M","-",turisme$Data)
gastos.ts<-ts(turisme[-1], start=c(2015,10), frequency = 12)

ef <- data.frame(x=1:86, y=turisme$`Eiv-Form`) #dades d'Eivissa i Formentera

serie_ef <- ts(ef$y,start=c(2015,10),frequency = 12)
plot.ts(serie_ef,  main="Eivissa i Formentera", xlab="Any", ylab="Despeses mensuals en €")

Les dades van des d’octubre de 2015 a novembre de 2022 (llavors tenim 7 cicles complets). Podem observar que l’efecte de la COVID és veu clarament en 2020.

Per veure millor l’estacionalitat de la sèrie visualitzem cada període mensualment:

seasonplot(serie_ef, col=c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),year.labels=TRUE, main="Estacionalitat d'Eivissa i Formentera", xlab="Mes", ylab=" Despeses en €")

legend("bottomright",
       legend = c(2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022),
       fill =  c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),      
       border = "black") 

Observam que efectivament, així com es comenta en el treball, el període de confinament va ser de principi de març, on comença a decréixer el valor dels preus, fins a finals de juny, on pareix que s’han tornat a recuperar els valors anteriors a la COVID.

Dividirem la sèrie en tres trams:

  • Des de 10-2015 fins 12-2018 (serie1_ef), amb el que predirem un període. És a dir predirem de 1-2019 a 12-2019. (amb aquest comprovam que els mètodes funcionen)

  • Des de 10-2015 fins 12-2019 (serie2_ef), amb la que predirem el període de la COVID (és a dir, l’any 2020)

  • Des de 10-2015 fins 12-2020 (serie3_ef), amb la que predirem el períodes després de la COVID (és a dir, l’any 2021)

Visualitzem-los:

serie1_ef <- ts(ef$y[1:39],start=c(2015,10),frequency = 12)
serie2_ef <- ts(ef$y[1:51],start=c(2015,10),frequency = 12)
serie3_ef <- ts(ef$y[1:63],start=c(2015,10),frequency = 12)

par(mfrow=c(2,2), mar=c(4,4,4,1)+.1)
plot.ts(serie_ef, main="Serie completa", xlab="Any", ylab="Euros")
plot.ts(serie1_ef, main="Serie abans de la COVID menys un cicle", xlab="Any", ylab="Euros")
plot.ts(serie2_ef, main="Serie abans de la COVID", xlab="Any", ylab="Euros")
plot.ts(serie3_ef, main="Serie abans i durant de la COVID", xlab="Any", ylab="Euros")


2. AJUST DELS MODELS

Abans d’aplicar algun model, estudiem un poc la sèrie:

serie_ef.lm <- lm(y~x, data=ef) 

plot.ts(ef$y, main = "Eivissa i Formentera", ylab = "Despeses mensuals en €", xlab="Índex de cada mes")

#dibuixam la recta de regressió sobre els punts
abline(serie_ef.lm, col='red')

Si dibuixam la recta de regressió sobre les nostres dades, tot i que aquestes estan molt disperses i no s’ajusten bé a la recta, podem observar que la tendència és una mica creixent, encara que és manté més o menys constant.

boxplot(serie_ef~cycle(serie_ef), xlab = "mesos", ylab = "Despeses en €", main="Boxplot d'Eivissa i Formentera")

Podem observar també la presència d’estacionalitat, que prenen els seus valors màxims a la temporada d’estiu i els seus mínims en l’hivern, fet que corresponen amb les dades turístiques a les illes.

Aplicarem diversos models per ajustar la nostra sèrie i fer-ne una predicció per llavors comparar quin és el millor.

Vegem com actua cada model:

2.1. ABANSE DE LA COVID


- REGRESSIÓ SEGMENTADA

Com hem comentat anteriorment la recta de regressió no s’ajusta bé a les dades, de fet el valor del \(R^2\) és molt baix:

serie1_ef.lm <- lm(y~x, data=ef[1:39,]) 
summary(serie1_ef.lm)$adj.r.squared
## [1] 0.02707681

Per això utilitzam el paquet segmented per ajustar a una recta de regressió a trossos.

Anem a utilitzar la funció selgmented() per veure quants de punts de canvi detecta:

punts_canvi_serie1_ef <-selgmented(serie1_ef.lm, Kmax = 10, type="bic")
## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 .. 
## 
## BIC to detect no. of breakpoints:
##        0        1        2        3        4        5        6        7 
## 419.0407 424.5861 429.3172 432.0721 430.7005 425.9556 386.0211 374.3229 
##        8     <NA>     <NA> 
## 373.3312 374.3312 375.3312 
## 
## No. of selected breakpoints: 6  (2 breakpoint(s) removed due to small slope diff)

Aplicam la funció segmented() que ens calcula la regressió segmentada:

serie1_ef.seg <- segmented(serie1_ef.lm, seg.Z=~x, psi = c(5,8,16,21,28,34))

Aquesta funció ens demana que introduïm els valors on es troben els punts de canvi, i aquesta ens dona el valor estimat. Aquests punts de canvi son:

summary(serie1_ef.seg)$psi
##        Initial      Est.    St.Err
## psi1.x       5  5.499227 0.4389605
## psi2.x       8 10.346005 0.4272700
## psi3.x      16 16.491929 0.4791051
## psi4.x      21 23.272538 0.4139897
## psi5.x      28 28.504180 0.3864705
## psi6.x      34 35.159760 0.3795130

Que corresponen, aproximadament, a: 2-2016, 7-2016, 1-2017, 8-2017, 2-2018 i 8-2018.

Obtenim que el valor de \(R^2\) per a la regressió segmentada és bastant alt

summary(serie1_ef.seg)$adj.r.squared
## [1] 0.8569661

Anem a visualitzar la regressió segmentada sobre les nostres dades

#per graficar-ho

p_serie1_ef <- ggplot(ef[1:39,], aes(x=x, y=y)) + geom_line()+
  ylim(c(0,1300)) +
  ggtitle("Regressió lineal i segmentada sobre les dades \nd'Eivissa i Formentera abans de la COVID") + 
  xlab('Índex del mes') + 
  ylab('Despeses mensuals en €')

my.coef1_ef <- coef(serie1_ef.lm) #coeficients de la recta de regressió lineal

p_serie1_ef <- p_serie1_ef + geom_abline(intercept=my.coef1_ef[1], slope=my.coef1_ef[2], color="green") #afegim la recta de regressió lineal

my.model1_ef <- data.frame(x=1:39, y=fitted(serie1_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada

p_serie1_ef <- p_serie1_ef + geom_line(data=my.model1_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines1_ef <- serie1_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats

p_serie1_ef <- p_serie1_ef + geom_vline(xintercept=my.lines1_ef, linetype="dashed") 

p_serie1_ef <- p_serie1_ef + theme(panel.background = element_blank())+ #eliminam el fons i la quadrícula del gràfic
  geom_vline(xintercept=0) + geom_hline(yintercept=0)

ggplotly(p_serie1_ef)

Visualment també es pot observar que la recta de regressió segmentada s’ajusta millor a les nostres dades.

Anem a calcular les equacions d’aquestes rectes, sabem que les rectes tenen la forma \(y=mx+n\), on \(m\) és la pendent i \(n\) el valor de tall de l’eix y.

Hi ha una funció del paquet segmented que ens dona aquesta valors de \(n\):

intercept(serie1_ef.seg)
## $x
##                Est.
## intercept1   933.81
## intercept2  -308.13
## intercept3  1969.50
## intercept4 -1013.20
## intercept5  4245.60
## intercept6 -2504.80
## intercept7  7206.40

També tenim una altra funció que ens calcula les pendents:

slope(serie1_ef.seg)
## $x
##            Est. St.Err. t value CI(95%).l CI(95%).u
## slope1  -92.588  24.405 -3.7938  -142.850   -42.325
## slope2  133.250  24.405  5.4599    82.988   183.510
## slope3  -86.896  18.449 -4.7102  -124.890   -48.901
## slope4   93.960  14.585  6.4423    63.922   124.000
## slope5 -132.000  24.405 -5.4088  -182.270   -81.741
## slope6  104.820  14.585  7.1867    74.779   134.860
## slope7 -171.390  34.514 -4.9656  -242.470  -100.300

Aleshores, la nostra recta segmentada és:

\(\hat{y}= \left\{ \begin{array}{lcc} -92.588x + 933.81 & si & x \leq \psi_1 \\ \\ 133.250x - 308.13 & si & \psi_1 < x \leq \psi_2 \\ \\ -86.896x + 1969.50 & si & \psi_2 < x \leq \psi_3 \\ \\ 93.960x - 1013.20 & si & \psi_3 < x \leq \psi_4 \\ \\ -132.000x + 4245.60 & si & \psi_4 < x \leq \psi_5 \\ \\ 104.820x - 2504.80 & si & \psi_5 < x \leq \psi_6 \\ \\ -171.390x + 7206.40 & si & \psi_6 < x \\ \end{array} \right.\)

on \(\psi_1= 5.5\), \(\psi_2= 10.35\), \(\psi_3= 16.49\), \(\psi_4=23.27\), \(\psi_5= 28.5\) i \(\psi_6= 35.16\).

Vegem els errors del model (ens interessa el RMSE):

accuracy(serie1_ef.seg)
##                        ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 4.372882e-15 61.79027 46.87676 -0.612519 6.148984 0.2692355

Anem a fer la predicció de 1-2019 a 12-2019.

El darrer punt de canvi que tenim és en agost de 2018 i, seguint el mateix criteri que amb les dades de Mallorca, podem considerar que els següents es donaran en 2-2019, 8-2019 i 2-2020. Necessitam calcular les pendents de les rectes d’entre agost de 2018 i febrer de 2020, per poder fer la predicció i calcular l’error.

Per calcular les pendents de les noves rectes farem la mitjana de les pendents anteriors. De les pendents ja calculades en el model obviarem la primera i la darrera, ja que no són vàlides. Per tant,

  • La pendent de 8-2018 a 2-2019 i de 8-2019 a 2-2020 serà : -109.448

  • La pendent de 2-2019 a 8-2019 serà : 110.677

Ara, seguint el mateix procediment que en el cas de Mallorca, els nous punts d’intersecció són: 5028.86, -3996.26 i 6349.61.

Llavors l’equació per a la predicció és:

\(\hat{y}= \left\{ \begin{array}{lcc} -109.448x + 5028.86 & si & x \leq \psi_7 \\ \\ 110.677x - 3996.26 & si & \psi_7 < x \leq \psi_8 \\ \\ -109.448x + 6349.61 & si & \psi_8 < x \\ \end{array} \right.\)

on \(\psi_7 = 41\) i \(\psi_8 = 47\).

# per graficar-la

p1_serie_ef <- ggplot(ef, aes(x=x, y=y)) + geom_line()+
  ylim(c(0,1300)) +
  ggtitle("Predicció d'e les Illes Balears'Eivissa i Formentera amb el \nmodel de regressió segmentada abans de la COVID") + 
  xlab('Índex del mes') + 
  ylab('Despeses mensuals en €')

my.model1_ef <- data.frame(x=1:39, y=fitted(serie1_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada

p1_serie_ef <- p1_serie_ef + geom_line(data=my.model1_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines1_ef <- serie1_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats

p1_serie_ef <- p1_serie_ef + geom_vline(xintercept=my.lines1_ef, linetype="dashed")

p1_serie_ef <- p1_serie_ef + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

p1_serie_ef <- p1_serie_ef + geom_abline(intercept = 5028.86, slope=-109.448, colour="green") +
                geom_abline(intercept = -3996.26, slope=110.677, colour="blue") +
                geom_abline(intercept = 6349.61, slope=-109.448, colour="orange")

ggplotly(p1_serie_ef)

Calculem l’error de la predicció:

o1_ef<-c(serie_ef[40:51]) #observacions reals de la sèrie de gener de 2019 a desembre de 2019

f0_ef <- -109.448 * 40 + 5028.68 #predicció de gener 2019

v1_ef=c(41:47)
f1_ef <- sapply(v1_ef, function(x) 110.677*x-3996.26) #predicció de febrer 2019 a agost 2019

v2_ef=c(48:51)
f2_ef <- sapply(v2_ef, function(x) -109.448*x + 6349.61 )#predicció de agost 2019 a desembre 2019

p1_ef<- c(f0_ef,f1_ef,f2_ef) #predicció de 2019

sqrt(sum((p1_ef-o1_ef)^2)/12) #error de la predicció
## [1] 109.3058


- DECOMPOSE

Recordem la sèrie

plot.ts(serie1_ef, main= "Eivissa i Formentera abans de la COVID", xlab="Any", ylab="Despeses mensuals en €")

Ja hem dit anteriorment que té una tendència mínima i podem observar també que no hi ha variabilitat.

El mètode clàssic de descomposició separa la sèrie en subseries corresponents a la tendència, la estacionalitat i el renou.

Aquestes components es poden combinar de manera additiva o multiplicativa. En el nostre cas utilitzam el model additiu: \(y_t = \mu_t + S_t + a_t\)

decompose_s1_ef<-decompose(serie1_ef)
plot(decompose_s1_ef, xlab="Any")

El decompose, per calcular aquestes noves tendències, el que fa és agafar les sis tendències anteriors i les sis posteriors de la sèrie original i en fa la mitjana. És per això que la primera que obtenim és en abril de 2016 i la darrera en juny de 2018. Notem que per a la predicció ens quedarem amb el valor de la darrera tendència del decompose.

t1_ef <- decompose_s1_ef$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t1_ef
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2015                                                                        
## 2016       NA       NA       NA 743.1446 747.8983 758.9400 759.6646 763.9542
## 2017 806.7158 817.0208 824.6025 837.5362 845.8621 843.2208 842.9546 837.8417
## 2018 831.0383 830.7542 836.4517 830.5658 816.6008 804.8000       NA       NA
##           Sep      Oct      Nov      Dec
## 2015                NA       NA       NA
## 2016 773.8896 782.0333 791.3829 797.4979
## 2017 833.9504 838.2292 840.2021 837.6517
## 2018       NA       NA       NA       NA

Els valors de les components estacionals els calcula fent la mitjana per mesos, és a dir, per calcular la componen de gener, agafa els valors de tots els geners que tenim a la sèrie original i en fa la mitjana. Per tant, només tenim 12 valors, un per a cada mes.

s1_ef <- decompose_s1_ef$seasonal
s1_ef <-  s1_ef[4:15] #estacionalitat de gener a desembre
s1_ef
##  [1] -306.83481 -201.66023 -246.34481  -95.47495   19.08352  117.80366
##  [7]  299.75769  236.97936  121.06727  197.11102  -22.20023 -119.28752

Anem a fer la predicció d’aquest model:

a1_ef <- c(s1_ef[7:12],s1_ef) #estacionalitat de juliol 2018 a desembre 2019 (es per poder fer la predicció)

pred1_decompose_ef <- sapply(a1_ef, function(x) 804.8000 + x) #predicció de juliol de 2018 a desembre 2019
p1_dec_ef<-c(serie1_ef[1:33], pred1_decompose_ef)

A1_ef<- data.frame("x" = ef[1:51,]$x, "y" = ef[1:51,]$y, "p"= p1_dec_ef)

grafica1_ef_dec <- ggplot(A1_ef)+
  geom_line(aes(x,p), color="red")+
  geom_line(aes(x,y))+
  geom_vline(xintercept=0) + geom_hline(yintercept=0)+  
  labs(title="Predicció d'Eivissa i Formntera abans de la COVID amb el model \nde descomposició", x="Índex del mes", y="Despeses mensuals en €")+ theme(panel.background = element_blank())

grafica1_ef_dec

Podem observar que aquest model fa una predicció bastant bona. Calculem l’error de la predicció:

(Notem que la predicció és des de juliol de 2018 a desembre de 2019 i per calcular l’error només volem de gener de 2019 a desembre de 2019.)

sqrt(sum((c(serie_ef[40:51]- pred1_decompose_ef[7:18]))^2)/12) #error predicció de gener a desembre de 2019
## [1] 84.1187


Notem que també tenim una altra instrucció a R per fer prediccions d’una sèrie, la funció predict(). Aquesta és basa en un model ETS, anem a veure la predicció:

prediccio1_ef <- predict(serie1_ef,h=12)
plot(prediccio1_ef, xlab="Any", ylab="Despeses mensuals en €")

Pareix que la predicció és bastant bona, ja que el cicle predit segueix un mateix patró que els anteriors. Anem a veure la predicció juntament amb la sèrie original:

df_prediccio1_ef <- data.frame(prediccio1_ef)

p1_ets_ef <- data.frame("x"= 1:51, "PointForecast"=c(serie1_ef,df_prediccio1_ef$Point.Forecast), "Lo80" = c(rep(NA,39),df_prediccio1_ef$Lo.80), "Hi80" = c(rep(NA,39),df_prediccio1_ef$Hi.80), "Lo95" = c(rep(NA,39),df_prediccio1_ef$Lo.95),"Hi95" = c(rep(NA,39),df_prediccio1_ef$Hi.95))

grafica_pred1_ets <- ggplot((ef[1:51,]))+
  geom_ribbon(data = p1_ets_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data =p1_ets_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = p1_ets_ef, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció del model ETS(M,N,M) a Eivissa i Formentera \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())

grafica_pred1_ets

Podem observar que la predicció és bastant bona, ja que continua seguint un mateix patró. I l’error comés és d’uns 76 euros.

accuracy(prediccio1_ef,serie2_ef)
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.9578372 60.39005 47.93061 -0.9764864 6.671517 0.4957655
## Test set     30.7451149 76.44673 66.73647  1.9682600 8.708841 0.6902821
##                   ACF1 Theil's U
## Training set 0.1630727        NA
## Test set     0.2101528 0.4449581


- SARIMA

Anem a veure quin model obtenim considerant un model estacional

Pel que hem vist anteriorment, podem considerar que no hi ha tendència, llavors no fa falta fer cap diferència a la part regular, no obstant, sí que cal fer una diferenciació d’orde 12. Vegem l’ACF i el PACF de la sèrie a predir:

par(mfrow=c(1,2))
acf(serie1_ef)
pacf(serie1_ef)

Per a la part regular obtenim un ARIMA(1,0,2) Feim una diferenciació a la part estacional, és a dir, d’ordre 12

serie1_ef_dif <- diff(serie1_ef,12)
plot(serie1_ef_dif, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")

Aquesta és la sèrie sense estacionalitat ni tendència, vegem com es modifiquen l’ACF i el PACF.

par(mfrow=c(2,1))
acf(serie1_ef_dif, lag=36)
pacf(serie1_ef_dif,lag=36)

Per a la part estacional obtenim que P=0, D=1, Q=0

És a dir, obtenim un ARIMA(1,0,2)(0,1,0)[12], vegem quin model els proposa R:

auto.arima(serie1_ef)
## Series: serie1_ef 
## ARIMA(0,1,1)(1,1,0)[12] 
## 
## Coefficients:
##           ma1     sar1
##       -0.6686  -0.4226
## s.e.   0.1826   0.2143
## 
## sigma^2 = 9792:  log likelihood = -156.79
## AIC=319.59   AICc=320.68   BIC=323.36

R ens proposa un ARIMA(0,1,1)(1,1,0)[12]. Per tant les propostes de model ARIMA són:

model1_ef<-arima(serie1_ef, order=c(1,0,2), seasonal = list(order=c(0,1,0), period=12))

model2_ef <- arima(serie1_ef, order=c(0,1,1), seasonal = list( order=c(1,1,0), period=12))

Amb uns errors d’uns 85 i 78 euros cada un.

accuracy(model1_ef)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -1.700455 84.64409 59.40385 -1.134031 7.905317 0.4795375
##                    ACF1
## Training set 0.04145267
accuracy(model2_ef)
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set -9.59489 77.62801 54.45449 -2.184592 7.553735 0.4395838 0.08070217

Anem a fer la predicció de 2019 d’aquests models:

fc_m1_ef <-forecast(model1_ef, h=12)
fc_m2_ef <- forecast(model2_ef,h=12)

Visualitzem-les

#Model 1
pre1_ef <- data.frame("Point Forecast" = serie1_ef, "Lo 80" = rep(NA,39), "Hi 80"=  rep(NA,39), "Lo 95" =  rep(NA,39), "Hi 95" =  rep(NA,39)) #dades abans de la predicció amb NA als intervals ja que només els volem per la predicció
 
pred_m1_ef <-data.frame(fc_m1_ef) 

grafica_m1_ef <- data.frame("x" = 1:51, "PointForecast" = c(pre1_ef$Point.Forecast,pred_m1_ef$Point.Forecast), "Lo80" = c(pre1_ef$Lo.80, pred_m1_ef$Lo.80), "Hi80"= c(pre1_ef$Hi.80, pred_m1_ef$Hi.80), "Lo95" = c(pre1_ef$Lo.95, pred_m1_ef$Lo.95), "Hi95" = c(pre1_ef$Hi.95, pred_m1_ef$Hi.95))

grafica1_ef <- ggplot(ef[1:51,])+
  geom_ribbon(data = grafica_m1_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m1_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m1_ef, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA(1,0,2)(0,1,0)[12] \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())

grafica1_ef

#Model 2

pred_m2_ef <-data.frame(fc_m2_ef)

grafica_m2_ef <- data.frame("x" = 1:51, "PointForecast" = c(pre1_ef$Point.Forecast,pred_m2_ef$Point.Forecast), "Lo80" = c(pre1_ef$Lo.80, pred_m2_ef$Lo.80), "Hi80"= c(pre1_ef$Hi.80, pred_m2_ef$Hi.80), "Lo95" = c(pre1_ef$Lo.95, pred_m2_ef$Lo.95), "Hi95" = c(pre1_ef$Hi.95, pred_m2_ef$Hi.95))


grafica2_ef <- ggplot(ef[1:51,])+
  geom_ribbon(data = grafica_m2_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m2_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m2_ef, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA (0,1,1)(1,1,0)[12] \nabans de la COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank())

grafica2_ef

Vegem i estudiem els errors de la predicció:

accuracy(fc_m1_ef, serie_ef[40:51], h=12)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -1.700455 84.64409 59.40385 -1.134031 7.905317 0.4795375
## Test set     52.755681 80.73996 64.25420  6.415168 8.130005 0.5186920
##                    ACF1
## Training set 0.04145267
## Test set             NA
accuracy(fc_m2_ef,serie_ef[40:51], h=12)
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set -9.59489 77.62801 54.45449 -2.184592 7.553735 0.4395838 0.08070217
## Test set     68.97582 93.75512 78.51529  8.037159 9.918771 0.6338147         NA
checkresiduals(fc_m1_ef)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,2)(0,1,0)[12]
## Q* = 5.3679, df = 5, p-value = 0.3727
## 
## Model df: 3.   Total lags used: 8
checkresiduals(fc_m2_ef)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(1,1,0)[12]
## Q* = 6.655, df = 6, p-value = 0.3539
## 
## Model df: 2.   Total lags used: 8
par(mfrow=c(1,2))
qqPlot(fc_m1_ef$residuals, id=FALSE, mean=mean(fc_m1_ef$residuals), sd=sd(fc_m1_ef$residuals))
qqPlot(fc_m2_ef$residuals, id=FALSE, mean=mean(fc_m2_ef$residuals), sd=sd(fc_m2_ef$residuals))

shapiro.test(fc_m1_ef$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m1_ef$residuals
## W = 0.93959, p-value = 0.03689
shapiro.test(fc_m2_ef$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m2_ef$residuals
## W = 0.94603, p-value = 0.06045
lillie.test(fc_m1_ef$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m1_ef$residuals
## D = 0.2282, p-value = 2.121e-05
lillie.test(fc_m2_ef$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m2_ef$residuals
## D = 0.19193, p-value = 0.0009128


Resum dels errors que cometen cada un dels models anteriors:

reg. segmentada descomposició clàssica ETS(M,N,M) ARIMA (1,0,2)(0,1,0)[12] ARIMA (0,1,1)(1,1,0)[12]
error model 61.79 60.39 84.64 77.62
error predicció 109.3058 84.1187 76.45 80.74 93.76


2.2. DURANT LA COVID


- REGRESSIÓ SEGMENTADA

El valor de \(R^2\) per a la regressió lineal és:

serie2_ef.lm <- lm(y~x, data= ef[1:51,])
summary(serie2_ef.lm)$adj.r.squared
## [1] 0.007558416

Aleshores utilitzem el paquet segmented per ajustar les nostres dades a una regressió lineal segmentada i millorar aquest valor. Vegem els punts de canvi:

punts_canvi_serie2_ef <-selgmented(serie2_ef.lm, Kmax = 10, type="bic")
## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 .. 
## 
## BIC to detect no. of breakpoints:
##        0        1        2        3        4        5        6        7 
## 551.5657 557.6175 564.7267 572.0761 573.9956 563.8713 564.6504 530.2193 
##        8        9       10 
## 508.0282 504.6689 508.9345 
## 
## No. of selected breakpoints: 8  (1 breakpoint(s) removed due to small slope diff)
serie2_ef.seg <- segmented(serie2_ef.lm, seg.Z=~x, psi = c(5,10,14,21,29,35,40,47))
summary(serie2_ef.seg)$psi
##        Initial      Est.    St.Err
## psi1.x       5  5.499227 0.4621615
## psi2.x      10 10.346005 0.4498531
## psi3.x      14 16.491929 0.5044278
## psi4.x      21 23.272538 0.4358708
## psi5.x      29 28.484784 0.4244957
## psi6.x      35 34.969280 0.4036198
## psi7.x      40 40.649312 0.3482145
## psi8.x      47 46.591245 0.3590303

Aquests són en: 3-2016, 7-2016, 1-2017, 8-2017, 1-2018, 8-2018, 2-2019 i 8-2019.

Ara el valor de \(R^2\) de la segmentada és

summary(serie2_ef.seg)$adj.r.squared
## [1] 0.852397

Anem a visualitzar la regressió segmentada sobre les nostres dades

p_serie2_ef <- ggplot(ef[1:51,], aes(x=x, y=y)) + geom_line()+
  ylim(c(0,1300))+
  ggtitle("Regressió lineal i segmentada sobre les dades \nd'Eivissa i Formentera durant la COVID") + xlab('índex del mes')+ ylab("Despeses mensuals en €")

my.coef2_ef <- coef(serie2_ef.lm) #coeficients de la recta de regressió lineal

p_serie2_ef <- p_serie2_ef + geom_abline(intercept=my.coef2_ef[1], slope=my.coef2_ef[2], color="green") #afegim la recta de regressió lineal

my.model2_ef <- data.frame(x=1:51, y=fitted(serie2_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada

p_serie2_ef <- p_serie2_ef + geom_line(data=my.model2_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines2_ef <- serie2_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats

p_serie2_ef <- p_serie2_ef + geom_vline(xintercept=my.lines2_ef, linetype="dashed") 

p_serie2_ef <- p_serie2_ef +  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

ggplotly(p_serie2_ef)

Per poder escriure la funció necessitam les pendents i els punts d’intersecció amb l’eix OY, que ens ho donen les següents funcions:

#punts d'intersecció
intercept(serie2_ef.seg)
## $x
##                Est.
## intercept1   933.81
## intercept2  -308.13
## intercept3  1969.50
## intercept4 -1013.20
## intercept5  4245.60
## intercept6 -2440.20
## intercept7  5850.90
## intercept8 -5229.00
## intercept9  7654.80
#pendents
slope(serie2_ef.seg)
## $x
##            Est. St.Err. t value CI(95%).l CI(95%).u
## slope1  -92.588  25.695 -3.6033  -144.870   -40.311
## slope2  133.250  25.695  5.1858    80.974   185.530
## slope3  -86.896  19.424 -4.4737  -126.410   -47.378
## slope4   93.960  15.356  6.1189    62.719   125.200
## slope5 -132.000  25.695 -5.1373  -184.280   -79.727
## slope6  102.710  19.424  5.2880    63.194   142.230
## slope7 -134.390  19.424 -6.9187  -173.900   -94.868
## slope8  138.190  19.424  7.1144    98.670   177.710
## slope9 -138.340  25.695 -5.3839  -190.620   -86.063

L’error del model de regressió segmentada és (ens interessa el RMSE):

accuracy(serie2_ef.seg)
##                         ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -3.342914e-15 65.36159 51.33711 -0.6744592 6.975136 0.2810544

Anem a fer la predicció, recordem que els punts de canvi són en 3-2016, 7-2016, 1-2017, 8-2017, 1-2018, 8-2018, 2-2019 i 8-2019. Aleshores, igual que abans els següents punts de canvi es donen en gener, agost i gener, per tant:

  • La pendent de 8-2019 a 1-2020 i de 8-2020 a 1-2021 serà: -117.762
  • La pendent de 1-2020 a 8-2020 serà : 117.0275

Seguin el mateix procediment que amb les dades de Mallorca, els nous punts d’intersecció són 6695.8, -5513.25 i 8339.33.

És a dir, la funció per a la predicció de 2020 és:

\(\hat{y}= \left\{ \begin{array}{lcc} -117.762x + 6695.8 & si & x \leq \psi_9 \\ \\ 117.0275x - 5513.25 & si & \psi_9 < x \leq \psi_{10} \\ \\ -117.762x + 8339.33& si & \psi_{10} < x \\ \end{array} \right.\)

on \(\psi_7 = 52\) i \(\psi_8 = 59\).

#ho graficam

p2_serie_ef <- ggplot(ef, aes(x=x, y=y)) + geom_line()+
  ylim(c(0,1400))+
  ggtitle("Predicció d'Eivissa i Formentera amb el model \nde regressió segmentada durant la COVID") + 
  xlab('Índex del mes') + 
  ylab('Despeses mensuals en €')

my.model2_ef <- data.frame(x=1:51, y=fitted(serie2_ef.seg)) #model nou amb els valors ajustats de la regressiño segmentada

p2_serie_ef <- p2_serie_ef + geom_line(data=my.model2_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines2_ef <- serie2_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats

p2_serie_ef <- p2_serie_ef + geom_vline(xintercept=my.lines2_ef, linetype="dashed") 

p2_serie_ef <- p2_serie_ef + theme(panel.background = element_blank()) + #eliminam el fons i la quadrícula del gràfic 
geom_vline(xintercept=0) + geom_hline(yintercept=0)

p2_serie_ef <- p2_serie_ef + geom_abline(intercept = 6695.8, slope=-117.762, colour="green") +
                geom_abline(intercept = -5513.25 , slope=117.0275, colour="blue") +
                geom_abline(intercept =8339.33 , slope=-117.762, colour="orange")

ggplotly(p2_serie_ef)

Calculem l’error de la predicció:

o2_ef<-c(serie_ef[52:63]) # dades reals per fer predicció del 2020

v3_ef=c(52:59)
f3_ef <- sapply(v3_ef, function(x) 117.0275*x-5513.25) #predicció de gener 2020 a agost 2020

v4_ef=c(60:63)
f4_ef <- sapply(v4_ef, function(x) -117.762*x + 8339.33) #predicció de setembre 2020 a desembre 2020

p2_ef <- c(f3_ef,f4_ef) #predicció de 2020

sqrt(sum((p2_ef-o2_ef)^2)/12)
## [1] 530.5246


- DECOMPOSE

La descomposició de la sèrie en aquest cas és la següent

decompose_s2_ef <- decompose(serie2_ef)
plot(decompose_s2_ef, xlab="Any")

On les components de tendència són:

t2_ef <- decompose_s2_ef$trend #tendències de la sèrie sense el tros a predir abans de la COVID 
t2_ef
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2015                                                                        
## 2016       NA       NA       NA 743.1446 747.8983 758.9400 759.6646 763.9542
## 2017 806.7158 817.0208 824.6025 837.5362 845.8621 843.2208 842.9546 837.8417
## 2018 831.0383 830.7542 836.4517 830.5658 816.6008 804.8000 801.1833 800.7592
## 2019 793.0312 797.2383 795.3000 794.5121 794.6762 798.6983       NA       NA
##           Sep      Oct      Nov      Dec
## 2015                NA       NA       NA
## 2016 773.8896 782.0333 791.3829 797.4979
## 2017 833.9504 838.2292 840.2021 837.6517
## 2018 794.4083 787.4583 785.2383 787.7225
## 2019       NA       NA       NA       NA

I les estacionals:

s2_ef <- decompose_s2_ef$seasonal #estacionalitat
s2_ef <-  s2_ef[4:15] #estacionalitat de gener a desembre

Fem la predicció:

a2_ef <- c(s2_ef[7:12],s2_ef) #estacionalitat de juliol 2019 a desembre 2029 (es per poder fer la predicció)

pred2_decompose_ef <- sapply(a2_ef, function(x) 798.6983 + x) #predicció de juliol de 2019 a desembre 2020

p2_dec_ef<-c(serie2_ef[1:45], pred2_decompose_ef)

A2_ef<- data.frame("x" = ef[1:63,]$x, "y" = ef[1:63,]$y, "p"= p2_dec_ef)

grafica2_ef_dec <- ggplot(A2_ef)+
  geom_line(aes(x,p), color="red")+
  geom_line(aes(x,y))+
  labs(title="Predicció d'Eivissa i Formentera durant la COVID amb el model \nde descomposició", x="Índex del mes", y="Despesese mensuals en €")+
  geom_vline(xintercept = 0)+ geom_hline(yintercept = 0)+ theme(panel.background = element_blank())

grafica2_ef_dec

L’error que es comet és:

sqrt(sum((c(serie_ef[52:63]- pred2_decompose_ef[7:18]))^2)/12)
## [1] 348.4049


Vegem, igual que abans, la predicció amb la funció predict():

prediccio2_ef <- predict(serie2_ef,h=12)
plot(prediccio2_ef, xlab="Any", ylab ="Despeses menusals en €")

Vegem com queda la predicció sobre la sèrie original:

df_prediccio2_ef <- data.frame(prediccio2_ef)
p2_ets_ef <- data.frame("x"= 1:63, "PointForecast"=c(serie2_ef,df_prediccio2_ef$Point.Forecast), "Lo80" = c(rep(NA,51),df_prediccio2_ef$Lo.80), "Hi80" = c(rep(NA,51),df_prediccio2_ef$Hi.80), "Lo95" = c(rep(NA,51),df_prediccio2_ef$Lo.95),"Hi95" = c(rep(NA,51),df_prediccio2_ef$Hi.95))

grafica_pred2_ets <- ggplot((ef[1:63,]))+
  geom_ribbon(data = p2_ets_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data =p2_ets_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = p2_ets_ef, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció del model ETS(A,N,A) a Eivissa i Formentera durant la COVID", y="Despeses mensuals en €", x="Índex del mes") + 
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0)+ theme(panel.background = element_blank())

grafica_pred2_ets

Si calculam l’error comés, aquest és d’uns 350 euros.

accuracy(prediccio2_ef, serie3_ef)
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set    1.253437  63.58661  49.70861 -0.6009612 6.749908 0.5855261
## Test set     -230.523743 350.31320 232.13630       -Inf      Inf 2.7343725
##                   ACF1 Theil's U
## Training set 0.1148920        NA
## Test set     0.4007914         0


- SARIMA

Quin model proposam nosaltres? Vegem l’ACF i el PACF:

par(mfrow=c(1,2))
acf(serie2_ef)
pacf(serie2_ef)

Per a la part regular obtenim que p=1, q=2 i d=0.

Observam que hi ha estacionalitat, llavors feim una diferenciació d’ordre 12.

serie2_ef_diff <- diff(serie2_ef,12)
plot(serie2_ef_diff, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")

Ara ja no s’observa l’estacionalitat, llavors hem de fer una diferenciació D=1. Vegem quins són els nous ACF i PACF.

par(mfrow=c(1,2))
acf(serie2_ef_diff,lag=36)
pacf(serie2_ef_diff,lag=36)

Obtenim que P=0 i Q=0. Per tant, el model que nosaltres proposam es un ARIMA(1,0,2)(0,1,0)[12]

R proposa el següent model:

auto.arima(serie2_ef)
## Series: serie2_ef 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.7560  -0.5981
## s.e.   0.1256   0.3089
## 
## sigma^2 = 7322:  log likelihood = -224.98
## AIC=455.95   AICc=456.66   BIC=460.87

Els dos models que tenim són

model3_ef <- arima(serie2_ef, order=c(1,0,2), seasonal = list(order=c(0,1,0), period=12))

model4_ef <- arima(serie2_ef, order=c(0,1,1), seasonal = list( order=c(0,1,1), period=12))

I els seus errors

accuracy(model3_ef)
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 6.692082 82.17198 60.48591 0.0511561 8.013875 0.4862884 0.02451003
accuracy(model4_ef)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -5.597112 71.89198 53.10005 -1.708372 7.395149 0.4269083
##                    ACF1
## Training set 0.08132473

Les prediccions són:

fc_m3_ef <- forecast(model3_ef, h=12)
fc_m4_ef <- forecast(model4_ef,h=12)
#primer model
pre2_ef <- data.frame("Point Forecast" = serie2_ef, "Lo 80" =  rep(NA,51), "Hi 80"= rep(NA,51), "Lo 95" = rep(NA,51), "Hi 95" = rep(NA,51))

pred_m3_ef <-data.frame(fc_m3_ef)

grafica_m3_ef <- data.frame("x" = 1:63, "PointForecast" = c(pre2_ef$Point.Forecast,pred_m3_ef$Point.Forecast), "Lo80" = c(pre2_ef$Lo.80, pred_m3_ef$Lo.80), "Hi80"= c(pre2_ef$Hi.80, pred_m3_ef$Hi.80), "Lo95" = c(pre2_ef$Lo.95, pred_m3_ef$Lo.95), "Hi95" = c(pre2_ef$Hi.95, pred_m3_ef$Hi.95))


grafica3_ef <- ggplot(ef[1:63,])+
  geom_ribbon(data = grafica_m3_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m3_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m3_ef, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Prediccó  d'Eivissa i Formentera amb el model ARIMA (1,0,2)(0,1,0)[12] \ndurant la COVID", x="Índex del mes", y="Despeses mensuals en €")+
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())


grafica3_ef

#segon model

pred_m4_ef <-data.frame(fc_m4_ef)

grafica_m4_ef <- data.frame("x" = 1:63, "PointForecast" = c(pre2_ef$Point.Forecast,pred_m4_ef$Point.Forecast), "Lo80" = c(pre2_ef$Lo.80, pred_m4_ef$Lo.80), "Hi80"= c(pre2_ef$Hi.80, pred_m4_ef$Hi.80), "Lo95" = c(pre2_ef$Lo.95, pred_m4_ef$Lo.95), "Hi95" = c(pre2_ef$Hi.95, pred_m4_ef$Hi.95))

grafica4_ef <- ggplot(ef[1:63,])+
  geom_ribbon(data = grafica_m4_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m4_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m4_ef, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA (0,1,1)(0,1,1)[12] \ndurant la COVID", x="Índex del mes", y="Despeses mensuals en €")+
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica4_ef

Vegem i estudiem els seus errors:

accuracy(fc_m3_ef,serie_ef[52:63], h=12)
##                       ME      RMSE       MAE       MPE     MAPE      MASE
## Training set    6.692082  82.17198  60.48591 0.0511561 8.013875 0.4862884
## Test set     -248.723052 374.77285 259.75035      -Inf      Inf 2.0883141
##                    ACF1
## Training set 0.02451003
## Test set             NA
accuracy(fc_m4_ef,serie_ef[52:63], h=12)
##                       ME      RMSE       MAE       MPE     MAPE      MASE
## Training set   -5.597112  71.89198  53.10005 -1.708372 7.395149 0.4269083
## Test set     -241.928797 360.06127 241.92880      -Inf      Inf 1.9450342
##                    ACF1
## Training set 0.08132473
## Test set             NA
checkresiduals(fc_m3_ef)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,2)(0,1,0)[12]
## Q* = 12.22, df = 7, p-value = 0.09355
## 
## Model df: 3.   Total lags used: 10
checkresiduals(fc_m4_ef)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
## Q* = 7.1127, df = 8, p-value = 0.5245
## 
## Model df: 2.   Total lags used: 10
par(mfrow=c(1,2))
qqPlot(fc_m3_ef$residuals, id=FALSE, mean=mean(fc_m3_ef$residuals), sd=sd(fc_m3_ef$residuals))
qqPlot(fc_m4_ef$residuals, id=FALSE, mean=mean(fc_m4_ef$residuals), sd=sd(fc_m4_ef$residuals))

shapiro.test(fc_m3_ef$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m3_ef$residuals
## W = 0.96971, p-value = 0.2151
shapiro.test(fc_m4_ef$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m4_ef$residuals
## W = 0.95558, p-value = 0.0541
lillie.test(fc_m3_ef$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m3_ef$residuals
## D = 0.15645, p-value = 0.003194
lillie.test(fc_m4_ef$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m4_ef$residuals
## D = 0.1427, p-value = 0.01114


Resum dels errors que comet cada model:

reg. segmentada descomposició clàssica ETS(A,N,A) ARIMA (1,0,2)(0,1,0)[12] ARIMA (0,1,1)(0,1,1)[12]
error model 65.36 63.59 82.17 71.84
error predicció 530.5246 348.405 350.31 374.77 360.06


2.3. DESPRÉS DE LA COVID


- REGRESSIÓ SEGMENTADA

El valor de \(R^2\) per a la regressió lineal és molt baix

serie3_ef.lm <- lm(y~x, data= ef[1:63,])
summary(serie3_ef.lm)$adj.r.squared
## [1] 0.00449053

Anem a fer ús del paquet segmented.

Vegem els punts de canvi:

punts_canvi_serie3_ef <-selgmented(serie3_ef.lm, Kmax = 10, type="bic")
## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 .. 
## 
## BIC to detect no. of breakpoints:
##        0        1        2        3        4        5        6        7 
## 698.6812 701.8040 707.7414 714.1851 721.0583 728.8687 713.0397 710.0976 
##        8        9       10 
## 714.0422 687.3246 694.5366 
## 
## No. of selected breakpoints: 8  (1 breakpoint(s) removed due to small slope diff)
serie3_ef.seg <- segmented(serie3_ef.lm, seg.Z=~x, psi = c(5,9,16,21,28,35,40,45,54,59))
summary(serie3_ef.seg)$psi
##         Initial      Est.    St.Err
## psi1.x        5  5.499227 0.5320942
## psi2.x        9 10.346005 0.5179006
## psi3.x       16 16.491929 0.5806260
## psi4.x       21 23.272538 0.5013678
## psi5.x       28 28.484784 0.4873200
## psi6.x       35 34.969280 0.4593709
## psi7.x       40 40.649312 0.3853447
## psi8.x       45 46.190155 0.3750704
## psi9.x       54 55.999956 0.2048643
## psi10.x      59 58.047110 0.2636145

Ens queden els punts de canvi a: 3-2016, 7-2016, 1-2017, 8-2017, 1-2018, 8-2018, 2-2019, 7-2019, 5-2020, 7-2020.

Ara, el valor de \(R^2\) de la segmentada és

summary(serie3_ef.seg)$adj.r.squared
## [1] 0.8531528

Anem a visualitzar la regressió segmentada damunt les nostres dades

p_serie3_ef <- ggplot(ef[1:63,], aes(x=x, y=y)) + geom_line()+
  ylim(c(0,1300))+
  ggtitle("Regressió lineal i segmentada sobre les dades  \nd'Eivissa i Formentera després de la COVID") + 
  xlab('índex del mes')+ 
  ylab('Despeses mensuals en €')

my.coef3_ef <- coef(serie3_ef.lm) #coeficients de la recta de regressió lineal

p_serie3_ef <- p_serie3_ef + geom_abline(intercept=my.coef3_ef[1], slope=my.coef3_ef[2], color="green") #afegim la recta de regressió lineal

my.model3_ef <- data.frame(x=1:63, y=fitted(serie3_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada

p_serie3_ef <- p_serie3_ef + geom_line(data=my.model3_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines3_ef <- serie3_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats

p_serie3_ef <- p_serie3_ef + geom_vline(xintercept=my.lines3_ef, linetype="dashed") 

p_serie3_ef <- p_serie3_ef + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

ggplotly(p_serie3_ef)

Per escriure la funció a trossos tenim:

#punts d'intersecció
intercept(serie3_ef.seg)
## $x
##                  Est.
## intercept1     933.81
## intercept2    -308.16
## intercept3    1969.60
## intercept4   -1013.80
## intercept5    4249.90
## intercept6   -2455.40
## intercept7    5932.00
## intercept8   -5595.60
## intercept9    6360.20
## intercept10 -23475.00
## intercept11   5317.20
#pendents
slope(serie3_ef.seg)
## $x
##             Est. St.Err. t value CI(95%).l CI(95%).u
## slope1   -92.589  29.584 -3.1297  -152.340   -32.843
## slope2   133.260  29.584  4.5043    73.509   193.000
## slope3   -86.907  22.363 -3.8861  -132.070   -41.743
## slope4    93.994  17.680  5.3165    58.289   129.700
## slope5  -132.180  29.584 -4.4680  -191.930   -72.436
## slope6   103.220  22.363  4.6155    58.055   148.380
## slope7  -136.630  22.363 -6.1097  -181.800   -91.469
## slope8   146.950  22.363  6.5712   101.790   192.120
## slope9  -111.880  12.078 -9.2638  -136.280   -87.494
## slope10  420.890  66.152  6.3625   287.290   554.480
## slope11  -75.126  29.584 -2.5394  -134.870   -15.380

L’error de la regressió segmentada és:

accuracy(serie3_ef.seg)
##                        ME    RMSE      MAE  MPE MAPE      MASE
## Training set 9.024968e-16 75.4706 55.09207 -Inf  Inf 0.2838393

Anem a fer la predicció per a l’any 2021. Recordem que els punts de canvi es donen a 3-2016, 7-2016, 1-2017, 8-2017, 1-2018, 8-2018, 2-2019, 7-2019, 5-2020, 7-2020. Degut a la pertorbació del COVID, sí que hi segueix havent un punt de canvi en estiu i l’altre en hivern successivament però ara no és donen al mateix mes. Per això, per predir els següents punts de canvi agafarem el mes més freqüent. Aleshores els següents punts de canvi seran en 1-2021, 7-2021 i 1-2022.

  • La pendent de 7-2020 a 1-2021 i de 7-2021 a 1-2022 serà : -116.899
  • La pendent de 1-2021 a 7-2021 serà : 179.663

Els nous punts d’intersecció es calculen de forma anàloga que a Mallorca. Per a les tres noves rectes obtenim que són \(n_1\)=7743.65, \(n_2\)=−11236.32 i \(n_3\)=9523.02.

Aleshores la predicció de 2021 serà:

\(\hat{y}= \left\{ \begin{array}{lcc} -116.899x + 7743.65 & si & x \leq \psi_{11} \\ \\ 179.663x - 11236.32 & si & \psi_{11} < x \leq \psi_{12} \\ \\ -116.899x + 9523.02 & si & \psi_{12} < x \\ \end{array} \right.\)

on \(\psi_7 = 64\) i \(\psi_8 = 70\).

Visualitzem-la:

p3_serie_ef <- ggplot(ef, aes(x=x, y=y)) + geom_line()+
  ylim(c(0,1300)) + 
  ggtitle("Predicció d'Eivissa i Formentera amb el model de \nregressió segmentada després de la COVID") + 
  xlab('índex del mes')+ 
  ylab('Despeses mensuals en €')

my.model3_ef <- data.frame(x=1:63, y=fitted(serie3_ef.seg)) #model nou amb els valors ajustats de la regressió segmentada

p3_serie_ef <- p3_serie_ef + geom_line(data=my.model3_ef, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines3_ef <- serie3_ef.seg$psi[,2]# guardam on estan els punts de canvi estimats

p3_serie_ef <- p3_serie_ef + geom_vline(xintercept=my.lines3_ef, linetype="dashed") 

p3_serie_ef <- p3_serie_ef + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

p3_serie_ef <- p3_serie_ef + geom_abline(intercept = 7743.65, slope=-116.899, colour="green") +
                geom_abline(intercept = -11236.32 , slope=179.663, colour="blue") +
                geom_abline(intercept = 9523.02 , slope=-116.899, colour="orange")

ggplotly(p3_serie_ef)

Vegem quin és aquest error que es comet:

o3_ef<-c(serie_ef[64:75]) # dades reals per fer predicció del 2021

v5_ef=c(64:70)
f5_ef <- sapply(v5_ef, function(x) 179.663*x-11236.32) #predicció de gener 2021 a juliol 2021

v6_ef=c(71:75)
f6_ef <- sapply(v6_ef, function(x) -116.899*x + 9523.02) #predicció d'agost 2021 a desembre 2021

p3_ef <- c(f5_ef,f6_ef) #predicció de 2021

sqrt(sum((p3_ef-o3_ef)^2)/13)
## [1] 192.499


- DECOMPOSE

Visualitzem la sèrie descomposada:

decompose_s3_ef <- decompose(serie3_ef)
plot(decompose_s3_ef, xlab="Any")

Les components de tendència són:

t3_ef <- decompose_s3_ef$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t3_ef
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2015                                                                        
## 2016       NA       NA       NA 743.1446 747.8983 758.9400 759.6646 763.9542
## 2017 806.7158 817.0208 824.6025 837.5362 845.8621 843.2208 842.9546 837.8417
## 2018 831.0383 830.7542 836.4517 830.5658 816.6008 804.8000 801.1833 800.7592
## 2019 793.0312 797.2383 795.3000 794.5121 794.6762 798.6983 802.9879 803.0183
## 2020 639.9054 615.4329 594.0237 577.0025 570.0400 569.5333       NA       NA
##           Sep      Oct      Nov      Dec
## 2015                NA       NA       NA
## 2016 773.8896 782.0333 791.3829 797.4979
## 2017 833.9504 838.2292 840.2021 837.6517
## 2018 794.4083 787.4583 785.2383 787.7225
## 2019 805.1829 776.5217 710.6350 662.4762
## 2020       NA       NA       NA       NA

I les d’estacionalitat:

s3_ef <- decompose_s3_ef$seasonal #estacionalitat
s3_ef <-  s3_ef[4:15] #estacionalitat de gener a desembre
s3_ef
##  [1] -226.71583 -201.30468 -200.48260 -185.77737  -83.92262  118.54038
##  [7]  314.53677  311.05855  156.86407  185.77125  -73.42270 -115.14520

Anem a fer la predicció:

a3_ef <- c(s3_ef[7:12],s3_ef) #estacionalitat de juliol 2020 a desembre 2021 (es per poder fer la predicció)

pred3_decompose_ef <- sapply(a3_ef, function(x) 569.5333  + x) #predicció de juliol de 2019 a desembre 2020

p3_dec_ef <-c(serie3_ef[1:57], pred3_decompose_ef)

A3_ef<- data.frame("x" = ef[1:75,]$x, "y" = ef[1:75,]$y, "p"= p3_dec_ef)

grafica3_ef_dec <- ggplot(A3_ef)+
  geom_line(aes(x,p), color="red")+
  geom_line(aes(x,y))+
labs(title="Predicció d'Eivissa i Formentera després de la COVID amb el model \nde descomposició", x="Índex del mes", y="Despeses mensualns en €")+
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica3_ef_dec

Calculem lerror de la predicció

sqrt(sum((c(serie_ef[64:75]- pred3_decompose_ef[7:18]))^2)/12)
## [1] 256.8545


Amb la funció predict() la predicció seria la següent:

prediccio3_ef <- predict(serie3_ef,h=12)
plot(prediccio3_ef, xlab="Any", ylab="Despeses mensuals en €")

Vegem-la juntament amb les nostres dades:

df_prediccio3_ef <- data.frame(prediccio3_ef)
p3_ets_ef <- data.frame("x"= 1:75, "PointForecast"=c(serie3_ef,df_prediccio3_ef$Point.Forecast), "Lo80" = c(rep(NA,63),df_prediccio3_ef$Lo.80), "Hi80" = c(rep(NA,63),df_prediccio3_ef$Hi.80), "Lo95" = c(rep(NA,63),df_prediccio3_ef$Lo.95),"Hi95" = c(rep(NA,63),df_prediccio3_ef$Hi.95))

grafica_pred3_ets <- ggplot((ef[1:75,]))+
  geom_ribbon(data = p3_ets_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data =p3_ets_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = p3_ets_ef, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció del model ETS(A,N,A) a Eivissa i Formentera \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+ 
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica_pred3_ets

Obtenim un error d’uns 138 euros.

accuracy(prediccio3_ef,serie_ef)
##                      ME     RMSE       MAE     MPE     MAPE      MASE      ACF1
## Training set  -2.165247 122.8461  85.26169    -Inf      Inf 0.6811543 0.1330175
## Test set     115.747684 137.8193 123.36967 15.6055 16.37964 0.9855984 0.5931860
##              Theil's U
## Training set        NA
## Test set      1.944982


- SARIMA

Quin model proposam nosaltres?Vegem l’ACF i el PACF:

par(mfrow=c(1,2))
acf(serie3_ef)
pacf(serie3_ef)

Per a la part regular: p=2, q=2 i d=0.

Sí hi ha indicació d’estacionalitat, llavors feim una diferenciació d’ordre 12 i tornam a calcular l’ACF i el PACF:

serie3_ef_diff <- diff(serie3_ef,12)
plot.ts(serie3_ef_diff, main="Sèrie sense estacionalitat", ylab="Sèrie diferenciada", xlab="Any")

par(mfrow=c(1,2))
acf(serie3_ef_diff)
pacf(serie3_ef_diff)

Aleshores obtenim un model ARIMA(2,0,2)(1,1,3)[12]

R proposa el següent model:

auto.arima(serie3_ef)
## Series: serie3_ef 
## ARIMA(0,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ma1      ma2     sma1
##       -0.1957  -0.5346  -0.7096
## s.e.   0.1395   0.1384   0.3900
## 
## sigma^2 = 20382:  log likelihood = -322
## AIC=651.99   AICc=652.88   BIC=659.64

Llavors tenim els següents models

model5_ef <- arima(serie3_ef, order=c(2,0,2), seasonal = list( order=c(1,1,3), period=12))

model6_ef <- arima(serie3_ef, order=c(0,1,2), seasonal = list( order=c(0,1,1), period=12))

Amb uns errors de:

accuracy(model5_ef)
##                     ME     RMSE      MAE  MPE MAPE      MASE       ACF1
## Training set -7.271777 109.7292 64.51072 -Inf  Inf 0.4889004 0.01064585
accuracy(model6_ef)
##                    ME    RMSE     MAE  MPE MAPE      MASE        ACF1
## Training set -13.9712 123.312 71.6025 -Inf  Inf 0.5426461 -0.03379418

I les seves prediccions són:

fc_m5_ef <- forecast(model5_ef, h=12)
fc_m6_ef <- forecast(model6_ef,h=12)
#primer model

pre3_ef <- data.frame("Point Forecast" = serie3_ef, "Lo 80" =  rep(NA,63), "Hi 80"= rep(NA,63), "Lo 95" = rep(NA,63), "Hi 95" = rep(NA,63))

pred_m5_ef <-data.frame(fc_m5_ef)


grafica_m5_ef <- data.frame("x" = 1:75, "PointForecast" = c(pre3_ef$Point.Forecast,pred_m5_ef$Point.Forecast), "Lo80" = c(pre3_ef$Lo.80, pred_m5_ef$Lo.80), "Hi80"= c(pre3_ef$Hi.80, pred_m5_ef$Hi.80), "Lo95" = c(pre3_ef$Lo.95, pred_m5_ef$Lo.95), "Hi95" = c(pre3_ef$Hi.95, pred_m5_ef$Hi.95))


grafica5_ef <- ggplot(ef[1:75,])+
  geom_ribbon(data = grafica_m5_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m5_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m5_ef, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA (2,0,2)(1,1,3)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+ 
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica5_ef

#segon model

pred_m6_ef <-data.frame(fc_m6_ef)

grafica_m6_ef <- data.frame("x" = 1:75, "PointForecast" = c(pre3_ef$Point.Forecast,pred_m6_ef$Point.Forecast), "Lo80" = c(pre3_ef$Lo.80, pred_m6_ef$Lo.80), "Hi80"= c(pre3_ef$Hi.80, pred_m6_ef$Hi.80), "Lo95" = c(pre3_ef$Lo.95, pred_m6_ef$Lo.95), "Hi95" = c(pre3_ef$Hi.95, pred_m6_ef$Hi.95))

grafica6_ef <- ggplot(ef[1:75,])+
  geom_ribbon(data = grafica_m6_ef, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m6_ef, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m6_ef, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció d'Eivissa i Formentera amb el model ARIMA (0,1,2)(0,1,1)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+ 
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica6_ef

Vegem quin és l’error de cada model i estudiem-los:

accuracy(fc_m5_ef,serie_ef[64:75], h=12)
##                      ME     RMSE       MAE      MPE     MAPE      MASE
## Training set  -7.271777 109.7292  64.51072     -Inf      Inf 0.4889004
## Test set     129.020379 164.9598 136.90003 17.47361 18.27666 1.0375094
##                    ACF1
## Training set 0.01064585
## Test set             NA
accuracy(fc_m6_ef,serie_ef[64:75], h=12)
##                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
## Training set -13.9712 123.3120  71.6025     -Inf      Inf 0.5426461 -0.03379418
## Test set     255.6337 270.3203 255.6337 33.14161 33.14161 1.9373431          NA
checkresiduals(fc_m5_ef)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,2)(1,1,3)[12]
## Q* = 3.0367, df = 5, p-value = 0.6943
## 
## Model df: 8.   Total lags used: 13
checkresiduals(fc_m6_ef)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,2)(0,1,1)[12]
## Q* = 5.8462, df = 10, p-value = 0.828
## 
## Model df: 3.   Total lags used: 13
par(mfrow=c(1,2))
qqPlot(fc_m5_ef$residuals, id=FALSE, mean=mean(fc_m5_ef$residuals), sd=sd(fc_m5_ef$residuals))
qqPlot(fc_m6_ef$residuals, id=FALSE, mean=mean(fc_m6_ef$residuals), sd=sd(fc_m6_ef$residuals))

shapiro.test(fc_m5_ef$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m5_ef$residuals
## W = 0.77159, p-value = 1.605e-08
shapiro.test(fc_m6_ef$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m6_ef$residuals
## W = 0.78424, p-value = 3.217e-08
lillie.test(fc_m5_ef$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m5_ef$residuals
## D = 0.15381, p-value = 0.0007865
lillie.test(fc_m6_ef$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m6_ef$residuals
## D = 0.16203, p-value = 0.0002898


Resum dels errors que comet cada un del models anteriors:

reg. segmentada descomposició clàssica ETS(A,N,A) ARIMA (2,0,2)(1,1,3)[12] ARIMA (0,1,2)(0,1,1)[12]
error model 75.47 122.85 109.73 123.31
error predicció 192.499 348.405 137.82 164.96 270.32